Victorian Forest Futures Population Vability Analysis: Central Highlands

1 Project setup

1.0.1 Packages

Install packages if necessary. Listed are installers for packages not on CRAN and the devlopment (most up-to-date) version of STEPS.

# library(devtools)
# install_github("steps-dev/steps")
# install_github("smwindecker/gdaltools")

Load packages

library(raster)
library(dplyr)
library(tibble)
library(sf)
library(rgdal)
library(readr)
library(readxl)
library(ggplot2)
library(lubridate)
library(magrittr)
library(tidyr)
library(foreach)
library(doMC)
library(future)
library(future.apply)
library(tidyr)
library(dismo)
library(gbm)
library(steps)
library(gdaltools)

1.0.2 Functions

source(file = "R/functions/pg.sf.R")
source(file = "R/functions/pg.pa.R")
source(file = "R/functions/read.vba.R")
source(file = "R/functions/proc.vba.R")
source(file = "R/functions/get.landis.vars.R")
source(file = "R/functions/rascc.R")
source(file = "R/functions/read.multi.line.header.R")
source(file = "R/functions/interpolate.climdat.R")
source(file = "R/functions/get.rst.prop.R")
source(file = "R/functions/get.rst.dat.R")

1.0.3 Paths

proj_path <- "/home/landis/rfst/"
# proj_path <- "D:/Users/ryan/Dropbox/Work/RFA_STEPS/rfst/"

2 Prepare variables

2.1 Analysis parameters

2.1.1 Timesteps

ntimesteps <- 50

ncores <- 20

2.2 Geographic data

2.2.1 Landscape

Set up base landscape layer

ch_rst <- raster(x = "data/grids/eco_v12.img") %T>%
  plot

ch_proj <- ch_rst@crs
ch_extent <- extent(ch_rst)
ch_res <- res(ch_rst)

2.2.2 Boundaries

rfa <- read_sf("data/shapefiles/RFA/")%>%
  st_transform(crs = ch_proj)
rfa

Simple feature collection with 81 features and 7 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -57985.97 ymin: 5665387 xmax: 763088.7 ymax: 6223446
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

pg.sf(rfa)

ch_rfa <- rfa[rfa$NAME== "CENTRAL HIGHLANDS",] 
ch_rfa

Simple feature collection with 1 feature and 7 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 318074.3 ymin: 5770575 xmax: 452175.7 ymax: 5906852
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

pg.sf(ch_rfa)

ch_mask <- rasterize(ch_rfa, ch_rst, field = 1, filename = "./output/landscape_vars/ch_mask.grd", overwrite = TRUE) 
ch_mask[ch_rst == 132] <- NA # removes lake areas from habitat
plot(ch_mask)

Think carefully about using this layer as projections are for zone 55, so distorts west of the state

victoria <- read_sf("data/shapefiles/vicstatepolygon/") %>%
   st_transform(crs = ch_proj)

2.3 Occurrence data

pseudo_abs <- read_sf("data/shapefiles/pseudo_absences/pseudo_absences.shp")%>%
  st_transform(crs = st_crs(ch_rst)) %>%
  mutate(PA = 0) %>%
  select(PA, geometry) %>%
  st_zm(drop = TRUE)
pseudo_abs

Simple feature collection with 35 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 325312 ymin: 5774360 xmax: 445518.4 ymax: 5901431
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

2.3.1 GG Occurrence data

Read in data

gg_ari <- read_excel(path = "data/tabular/BoA_SB_Combined_VBA_upload_v4.xls") %>%
  dplyr::select(-starts_with("leave")) %>%
  rename("lon" = `X-coordinate (easting or longitude)`, "lat" = `Y-coordinate (northing or latitude)`) %>%
   fill(lon, lat, .direction = "down") %>%
  filter(`Taxon Name` == "Misc Target taxa not found") %>%
  select(lon, lat) %>%
  mutate(PA = 0) %>%
  st_as_sf(coords = c("lon", "lat"), crs = st_crs(28355)) %>%
  st_transform(crs = st_crs(ch_rst))

-
/
                                                                                                                  

New names:
* `leave blank` -> `leave blank...2`
* `leave blank` -> `leave blank...3`
* `leave blank` -> `leave blank...6`
* `leave blank` -> `leave blank...7`
* `leave blank` -> `leave blank...25`

gg_ari

Simple feature collection with 51 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 333214 ymin: 5779763 xmax: 451285 ymax: 5937215
epsg (SRID):    28355
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

gg_vba <- proc.vba("data/tabular/vba_gg_all_20190509.csv", project.crs = ch_proj)
gg_vba

Simple feature collection with 3086 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 248825.6 ymin: 5725915 xmax: 742929.5 ymax: 6008846
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

gg_pa_all <- gg_vba %>%
  rbind(gg_ari)
gg_pa_all

Simple feature collection with 3137 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 248825.6 ymin: 5725915 xmax: 742929.5 ymax: 6008846
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

plot_gg_pa_all <- pg.pa(gg_pa_all, ch_rfa) +
  geom_sf(data = victoria, fill = NA)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

plot_gg_pa_all

gg_pa_ch <- gg_pa_all[ch_rfa,]
gg_pa_ch

Simple feature collection with 1308 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 333053.7 ymin: 5783845 xmax: 448916.2 ymax: 5887474
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

gg_pa_ch_psab <- gg_pa_ch %>%
  rbind(pseudo_abs)

st_write(gg_pa_ch, "output/shapefiles/pa/ch/gg_pa_ch.shp", delete_layer = TRUE)

Deleting layer `gg_pa_ch' using driver `ESRI Shapefile'
Writing layer `gg_pa_ch' to data source `output/shapefiles/pa/ch/gg_pa_ch.shp' using driver `ESRI Shapefile'
features:       1308
fields:         1
geometry type:  Point

st_write(gg_pa_ch_psab, "output/shapefiles/pa/ch/gg_pa_ch_psab.shp", delete_layer = TRUE)

Deleting layer `gg_pa_ch_psab' using driver `ESRI Shapefile'
Writing layer `gg_pa_ch_psab' to data source `output/shapefiles/pa/ch/gg_pa_ch_psab.shp' using driver `ESRI Shapefile'
features:       1343
fields:         1
geometry type:  Point

pg_gg_pa_ch <- pg.pa(gg_pa_ch %>% mutate(pseudo_abs = "vba") %>% rbind(pseudo_abs %>% mutate(pseudo_abs = "pseudo_absences")), ch_rfa) +
  facet_grid(.~pseudo_abs)
pg_gg_pa_ch

2.3.2 LBP Occurrence data

lb_pa_all <- proc.vba("data/tabular/vba_lb_all_20190509.csv", project.crs = ch_proj)

plot_lb_pa_all <- pg.pa(lb_pa_all, ch_rfa) +
  geom_sf(data = victoria, fill = NA)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

plot_lb_pa_all

lb_pa_ch <- lb_pa_all[ch_rfa,]
lb_pa_ch

Simple feature collection with 1404 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 364618.3 ymin: 5798912 xmax: 447312.4 ymax: 5868055
epsg (SRID):    NA
proj4string:    +proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs

lb_pa_ch_psab <- lb_pa_ch %>%
  rbind(pseudo_abs)

st_write(lb_pa_ch, "output/shapefiles/pa/ch/lb_pa_ch.shp", delete_layer = TRUE)

Deleting layer `lb_pa_ch' using driver `ESRI Shapefile'
Writing layer `lb_pa_ch' to data source `output/shapefiles/pa/ch/lb_pa_ch.shp' using driver `ESRI Shapefile'
features:       1404
fields:         1
geometry type:  Point

st_write(lb_pa_ch_psab, "output/shapefiles/pa/ch/lb_pa_ch_psab.shp", delete_layer = TRUE)

Deleting layer `lb_pa_ch_psab' using driver `ESRI Shapefile'
Writing layer `lb_pa_ch_psab' to data source `output/shapefiles/pa/ch/lb_pa_ch_psab.shp' using driver `ESRI Shapefile'
features:       1439
fields:         1
geometry type:  Point

pg_lb_pa_ch <- pg.pa(lb_pa_ch %>% mutate(pseudo_abs = "vba") %>% rbind(pseudo_abs %>% mutate(pseudo_abs = "pseudo_absences")), ch_rfa) +
  facet_grid(.~pseudo_abs)
pg_lb_pa_ch

save.image("output/session_images/landscape-occupancy.RData")

#load(file = "output/session_images/landscape-occupancy.RData")

2.4 Landis output

2.4.1 Scenario 1: Current planned burning, current timber harvesting, climate change under RCP 4.5, current bushfire regime (Buisness as usual)

aa <- get.landis.vars(
  scn_path = "~/landis_ch_s1_pb-th-cc_50/",
  proj_path = proj_path,
  out_path = "/output/habitat_vars/",
  scn_id = "s1",
  proj_mask = ch_mask,
  timesteps = 3,
  cores = 4
  )

plot(s1_vars_landis[[ntimesteps + 1]])

save(s1_vars_landis, file = "output/habitat_vars/s1_vars_landis")

# load(file = "output/habitat_vars/s1_vars_landis_fut")

2.4.2 Scenario 2: Current planned burning, current timber harvesting, climate change under RCP 8.5, increatred bushfire regime

2.4.3 Scenario 3: Current planned burning, current timber harvesting, no climate change, current bushfire regime

s3_vars_landis <- get.landis.vars(
  scn_path = "~/landis_ch_s3_pb-th-cc0_50/",
  proj_path = proj_path,
  out_path = "/output/habitat_vars/",
  scn_id = "s3",
  proj_mask = ch_mask,
  timesteps = ntimesteps,
  cores = ncores
  )

plot(s3_vars_landis[[1]])

plot(s3_vars_landis[[ntimesteps + 1]])

save(s3_vars_landis, file = "output/habitat_vars/s3_vars_landis")
# load(file = "output/habitat_vars/s3_vars_landis_fut")

2.4.4 Scenario 4: Current planned burning, no harvesting, climate change under RCP 4.5, current bushfire regime

s4_vars_landis <- get.landis.vars(
  scn_path = "~/landis_ch_s4_pb-th0-cc_50/",
  proj_path = proj_path,
  out_path = "/output/habitat_vars/",
  scn_id = "s4",
  proj_mask = ch_mask,
  timesteps = ntimesteps,
  cores = ncores,
  harvest_timber = FALSE
  )

plot(s4_vars_landis[[1]])

plot(s4_vars_landis[[ntimesteps + 1]])

save(s4_vars_landis, file = "output/habitat_vars/s4_vars_landis")
# load(file = "output/habitat_vars/s4_vars_landis_fut")

2.5 ARI Data

2.5.1 Anisotronic heating x ruggedness

ahr <-raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/Anisotrophic_Heating_Ruggedness") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_ahr.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.2 Log verticical distance all wetlands

lvdaw <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_all_wetlands_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_lvdaw.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.3 Log vertical distance major streams

lvdma <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_major_streams_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_lvdma.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.4 Log vertical distance minor streams

lvdmi <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_minor_streams_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask) %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.5 Log vertical distance saline wetlands

lvdsw <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/log_vertical_distance_saline_wetlands_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_lvdsw.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.6 Relative annual insolation

rai <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/relative_annual_insolation_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_rai.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.7 Thorium

tho <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/sept2014thorium") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_tho.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.8 Visible sky index

vsky <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/visible_sky_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_vsky.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.9 Topographic wetness index

twi <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/wetness_index_topocrop_sept2012") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/habitat_vars/ari_twi.grd") %T>%
  plot

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.5.10 Wind exposition

Problem with layer

# winex <- raster(x = "data/grids/Env_covariates_noClimate_VicGrid94/wind_exposition2015") %>%
  # projectRaster(to = ch_mask) %>%
  # mask(mask = ch_mask, filename = "output/habitat_vars/ari_winex.grd") %T>%
  # plot

2.5.11 ARI variables stack

2.5.11.1 Initial variables

ari_iv <- stack(lvdaw, lvdma, lvdmi, lvdsw)
names(ari_iv) <- c("lvdaw", "lvdma", "lvdmi", "lvdsw")

2.5.11.2 Future variables

Repeates into list for each year

ari_v <- vector("list", ntimesteps + 1)
for(i in 1:(ntimesteps+1)){
  ari_v[[i]] <- ari_iv
}

2.6 Climactic data

pc ~ percentage change ac ~ absolute change

tmax ~ max temperature tmin ~ minimum temperature prec ~ precipitation evtr ~ evapotranspiration

01 ~ January 07 ~ July

4.5 ~ rcp 4.5 8.5 ~ rcp 8.5

2.6.1 Current climate

From WorldClim

2.6.1.1 Precipitation in January

prec01 <- raster("data/grids/wc2.0_30s_prec_01.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/prec01.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(prec01) <- "prec01"
plot(prec01)

2.6.1.2 Precipitation in July

prec07 <- raster("data/grids/wc2.0_30s_prec_07.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/prec07.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(prec07) <- "prec07"
plot(prec07)

2.6.1.3 Max temperature in January

tmax01 <- raster("data/grids/wc2.0_30s_tmax_01.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/tmax01.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(tmax01) <- "tmax01"
plot(tmax01)

2.6.1.4 Minimum temperature in July

tmin07 <- raster("data/grids/wc2.0_30s_tmin_07.tif") %>%
  projectRaster(to = ch_mask) %>%
  mask(mask = ch_mask, filename = "output/clim_vars/tmin07.grd")

3 projected point(s) not finiteno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

names(tmin07) <- "tmin07"
plot(tmin07)

2.6.2 Future climate

Data from Climate Change in Australia

2.6.2.1 Change data

There should be new data to replace this ##### Absolute change in temperature Get raw data

#absolute change in temp
raw_tmax01_4.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmax_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmax_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

raw_tmax01_8.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmax_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmax_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

raw_tmin07_4.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmin_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmin_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]") 

Read 4 items

raw_tmin07_8.5_ac <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/tasmin_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_abs-change-wrt-seasavg-clim_native.csv?tasmin_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

Years data available

n_tmax01_4.5 <- as.numeric(sub("-.*", "", unique(raw_tmax01_4.5_ac$time)))
n_tmax01_8.5 <- as.numeric(sub("-.*", "", unique(raw_tmax01_8.5_ac$time)))
n_tmin07_4.5 <- as.numeric(sub("-.*", "", unique(raw_tmin07_4.5_ac$time)))
n_tmin07_8.5 <- as.numeric(sub("-.*", "", unique(raw_tmin07_8.5_ac$time)))

Reprojected layers

tmax01_4.5_ac <- rascc(raw_tmax01_4.5_ac, new.proj.layer = ch_mask, filename = "output/clim_vars/tmax01_4.5_ac.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

tmax01_8.5_ac <- rascc(raw_tmax01_8.5_ac, new.proj.layer = ch_mask, filename = "output/clim_vars/tmax01_8.5_ac.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

tmin07_4.5_ac <- rascc(raw_tmin07_4.5_ac, new.proj.layer = ch_mask, filename = "output/clim_vars/tmin07_4.5_ac.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

tmin07_8.5_ac <- rascc(raw_tmin07_8.5_ac, new.proj.layer = ch_mask, filename = "output/clim_vars/tmin07_8.5_ac.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.6.2.1.1 Percentage change in precipitation

Get raw data

raw_prec01_4.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

raw_prec01_8.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_january[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

raw_prec07_4.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp45_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

raw_prec07_8.5_pc <- read.multi.line.header(file = "http://nrm-erddap.nci.org.au/erddap/griddap/pr_Amon_CSIRO-Mk3-6-0_rcp85_r1i1p1_perc-change-wrt-seassum-clim_native.csv?pr_july[(2025-01-01T12:00:00Z):1:(2090-01-01T12:00:00Z)][(-40.10297775):1:(-32.64199448)][(140.625):1:(150)]")

Read 4 items

Years data available

n_prec01_4.5 <- as.numeric(sub("-.*", "", unique(raw_prec01_4.5_pc$time)))
n_prec01_8.5 <- as.numeric(sub("-.*", "", unique(raw_prec01_4.5_pc$time)))
n_prec07_4.5 <- as.numeric(sub("-.*", "", unique(raw_prec07_4.5_pc$time)))
n_prec07_8.5 <- as.numeric(sub("-.*", "", unique(raw_prec07_4.5_pc$time)))

Reprojected layers

prec01_4.5_pc <- rascc(raw_prec01_4.5_pc, new.proj.layer = ch_mask, filename = "output/clim_vars/prec01_4.5_pc.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

prec01_8.5_pc <- rascc(raw_prec01_8.5_pc, new.proj.layer = ch_mask, filename = "output/clim_vars/prec01_8.5_pc.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

prec07_4.5_pc <- rascc(raw_prec07_4.5_pc, new.proj.layer = ch_mask, filename = "output/clim_vars/prec07_4.5_pc.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

prec07_8.5_pc <- rascc(raw_prec07_8.5_pc, new.proj.layer = ch_mask, filename = "output/clim_vars/prec07_8.5_pc.grd")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

2.6.2.2 Absolute predicted values

2.6.2.2.1 Temperature

add new function that writes these climate layers to files.

2.6.2.2.1.1 Jan max temperature RCP 4.5

tmax01_4.5 <- rst.op(input1 = tmax01,
                     input2 = tmax01_4.5_ac,
                     proj_mask = ch_mask,
                     op = "addabs",
                     filename = "output/clim_vars/tmax01_4.5",
                     layernames = n_tmax01_4.5) %T>%
  plot

2.6.2.2.1.2 Jan max temperature RCP 8.5

tmax01_8.5 <- rst.op(input1 = tmax01,
                     input2 = tmax01_8.5_ac,
                     proj_mask = ch_mask,
                     op = "addabs",
                     filename = "output/clim_vars/tmax01_8.5",
                     layernames = n_tmax01_8.5) %T>%
  plot

2.6.2.2.1.3 July min temperature RCP 4.5

tmin07_4.5 <- rst.op(input1 = tmin07,
                     input2 = tmin07_4.5_ac,
                     proj_mask = ch_mask,
                     op = "addabs",
                     filename = "output/clim_vars/tmin07_4.5",
                     layernames = n_tmin07_4.5) %T>%
  plot

2.6.2.2.1.4 July min temperature RCP 8.5

tmin07_8.5 <- rst.op(input1 = tmin07,
                     input2 = tmin07_8.5_ac,
                     proj_mask = ch_mask,
                     op = "addabs",
                     filename = "output/clim_vars/tmin07_8.5",
                     layernames = n_tmin07_8.5) %T>%
  plot

2.6.2.2.2 Precipitation

2.6.2.2.3 January precipitation RCP 4.5

prec01_4.5 <- rst.op(input1 = prec01,
                     input2 = prec01_4.5_pc,
                     proj_mask = ch_mask,
                     op = "addper",
                     filename = "output/clim_vars/prec01_4.5",
                     layernames = n_prec01_4.5) %T>%
  plot

2.6.2.2.4 January precipitation RCP 8.5

prec01_8.5 <- rst.op(input1 = prec01,
                     input2 = prec01_8.5_pc,
                     proj_mask = ch_mask,
                     op = "addper",
                     filename = "output/clim_vars/prec01_8.5",
                     layernames = n_prec01_8.5) %T>%
  plot

2.6.2.2.5 July precipitation RCP 4.5

prec07_4.5 <- rst.op(input1 = prec07,
                     input2 = prec07_4.5_pc,
                     proj_mask = ch_mask,
                     op = "addper",
                     filename = "output/clim_vars/prec07_4.5",
                     layernames = n_prec07_4.5) %T>%
  plot

2.6.2.2.6 July precipitation RCP 8.5

prec07_8.5 <- rst.op(input1 = prec07,
                     input2 = prec07_8.5_pc,
                     proj_mask = ch_mask,
                     op = "addper",
                     filename = "output/clim_vars/prec07_8.5",
                     layernames = n_prec07_8.5) %T>%
  plot

2.6.2.3 Interploate prediction data

Currently just repeats predictions for nearest year, need updating to weight predictions within 10 years, cf data from CSIRO

2.6.3 Climate variable sets

2.6.3.1 Initial climate variables

clim_iv <- stack(prec01, prec07, tmax01, tmin07)

2.6.3.2 Future climate variables

clim_fv_cc0 <- vector("list", 50)
for(i in 1:50){
  clim_fv_cc0[[i]] <- clim_iv
}


clim_fv_4.5 <- mapply(prec01_4.5_int, prec07_4.5_int, tmax01_4.5_int, tmin07_4.5_int, FUN = stack)
clim_fv_8.5 <- mapply(prec01_8.5_int, prec07_8.5_int, tmax01_8.5_int, tmin07_8.5_int, FUN = stack)

save.image(file = "output/input_variables.RData")

2.7 Distribution model variables

2.7.1 Species LANDIS variables

2.7.1.1 GG LANDIS vars

gglv <- c("ggf", "ggd", "hbt_1k")

2.7.1.2 LB LANDIS vars

lblv <- c("lbm", "hbt_3h")

2.7.2 Scenario 1

2.7.2.1 Greater Glider

2.7.2.1.1 Scenario 1 initial variables

s1_iv_gg <- stack(s1_vars_landis_init[[c(gglv)]], clim_iv)

save(s1_iv_gg, file = "output/habitat_vars/s1_iv_gg")

2.7.2.1.2 Scenaio 1 model data

s1_md_gg <- s1_iv_gg %>%
  raster::extract(y = st_coordinates(gg_pa_ch)) %>%
  cbind("PA" = gg_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s1_md_gg

save(s1_md_gg, file = "output/habitat_vars/s1_md_gg")

2.7.2.1.3 Scenario 1 future variables

s1_fv_gg <- s1_vars_landis_fut %>%
  sapply("[[", gglv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s1_fv_gg, file = "output/habitat_vars/s1_fv_gg")

2.7.2.2 Leadbeater’s Possum

2.7.2.2.1 Scenario 1 initial variables

s1_iv_lb <- stack(s1_vars_landis_init[[c(lblv)]], clim_iv)

save(s1_iv_lb, file = "output/habitat_vars/s1_iv_lb")

2.7.2.2.2 Scenaio 1 model data

s1_md_lb <- s1_iv_lb %>%
  raster::extract(y = st_coordinates(lb_pa_ch)) %>%
  cbind("PA" = lb_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s1_md_lb

save(s1_md_lb, file = "output/habitat_vars/s1_md_lb")

2.7.2.2.3 Scenario 1 future variables

s1_fv_lb <- s1_vars_landis_fut %>%
  sapply("[[", lblv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s1_fv_lb, file = "output/habitat_vars/s1_fv_lb")

2.7.3 Scenario 2

awaiting updated climate variables

2.7.4 Scenario 3

2.7.4.1 Greater Glider

2.7.4.1.1 Scenario 3 initial variables

s3_iv_gg <- stack(s3_vars_landis_init[[c(gglv)]], clim_iv)

save(s3_iv_gg, file = "output/habitat_vars/s3_iv_gg")

2.7.4.1.2 Scenaio 3 model data

s3_md_gg <- s3_iv_gg %>%
  raster::extract(y = st_coordinates(gg_pa_ch)) %>%
  cbind("PA" = gg_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s3_md_gg

save(s3_md_gg, file = "output/habitat_vars/s3_md_gg")

2.7.4.1.3 Scenario 3 future variables

s3_fv_gg <- s3_vars_landis_fut %>%
  sapply("[[", gglv) %>%
  mapply(clim_fv_cc0, FUN = stack)

save(s3_fv_gg, file = "output/habitat_vars/s3_fv_gg")

2.7.4.2 Leadbeater’s Possum

2.7.4.2.1 Scenario 3 initial variables

s3_iv_lb <- stack(s3_vars_landis_init[[c(lblv)]], clim_iv)

save(s3_iv_lb, file = "output/habitat_vars/s3_iv_lb")

2.7.4.2.2 Scenaio 3 model data

s3_md_lb <- s3_iv_lb %>%
  raster::extract(y = st_coordinates(lb_pa_ch)) %>%
  cbind("PA" = lb_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s3_md_lb

save(s3_md_lb, file = "output/habitat_vars/s3_md_lb")

2.7.4.2.3 Scenario 3 future variables

s3_fv_lb <- s3_vars_landis_fut %>%
  sapply("[[", lblv) %>%
  mapply(clim_fv_cc0, FUN = stack)

save(s3_fv_lb, file = "output/habitat_vars/s3_fv_lb")

2.7.5 Scenario 4

2.7.5.1 Greater Glider

2.7.5.1.1 Scenario 4 initial variables

s4_iv_gg <- stack(s4_vars_landis_init[[c(gglv)]], clim_iv)

save(s4_iv_gg, file = "output/habitat_vars/s4_iv_gg")

2.7.5.1.2 Scenaio 4 model data

s4_md_gg <- s4_iv_gg %>%
  raster::extract(y = st_coordinates(gg_pa_ch)) %>%
  cbind("PA" = gg_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s4_md_gg

save(s4_md_gg, file = "output/habitat_vars/s4_md_gg")

2.7.5.1.3 Scenario 4 future variables

s4_fv_gg <- s4_vars_landis_fut %>%
  sapply("[[", gglv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s4_fv_gg, file = "output/habitat_vars/s4_fv_gg")

2.7.5.2 Leadbeater’s Possum

2.7.5.2.1 Scenario 4 initial variables

s4_iv_lb <- stack(s4_vars_landis_init[[c(lblv)]], clim_iv)

save(s4_iv_lb, file = "output/habitat_vars/s4_iv_lb")

2.7.5.2.2 Scenaio 4 model data

s4_md_lb <- s4_iv_lb %>%
  raster::extract(y = st_coordinates(lb_pa_ch)) %>%
  cbind("PA" = lb_pa_ch$PA, .) %>%
  as.data.frame %>%
  na.omit 

s4_md_lb

save(s4_md_lb, file = "output/habitat_vars/s4_md_lb")

2.7.5.2.3 Scenario 4 future variables

s4_fv_lb <- s4_vars_landis_fut %>%
  sapply("[[", lblv) %>%
  mapply(clim_fv_4.5, FUN = stack)

save(s4_fv_lb, file = "output/habitat_vars/s4_fv_lb")

3 Distribution model fit and prediction

3.1 Greater Glider

3.1.1 Scenario 1

s1_mod_gg <- gbm.step(data = s1_md_gg, gbm.x = 2:8, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s1_mod_gg)

#s1_ip_gg <- predict(s1_iv_gg, s1_mod_gg, type = "response", n.trees = s1_mod_gg$gbm.call$best.trees, filename = "output/habitat_pred/s1_fp_gg-000.tif")

#writeRaster(s1_ip_gg, "output/habitat_pred/s1_ip_gg.tif", overwrite = TRUE)
s1_ip_gg <- raster(x = "output/habitat_pred/s1_ip_gg.tif")

plot(s1_ip_gg)

registerDoMC(cores = 20)

s1_fp_gg <- foreach(i = seq_len(length(s1_fv_gg))) %dopar% {
  predict(s1_fv_gg[[i]],
          s1_mod_gg,
          type = "response",
          n.trees = s1_mod_gg$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s1_fp_gg-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s1_fp_gg <- stack(s1_fp_gg)

save(s1_fp_gg, file = "output/habitat_pred/s1_fp_gg")
writeRaster(s1_fp_gg, filename = "output/habitat_pred/s1_fp_gg.tif")

3.1.2 Scenario 3

s3_mod_gg <- gbm.step(data = s3_md_gg, gbm.x = 3:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s3_mod_gg)

 # s3_ip_gg <- predict(s3_iv_gg, s3_mod_gg, type = "response", n.trees = s3_mod_gg$gbm.call$best.trees, filename ="output/habitat_pred/s3_fp_gg-000.tif")

 # writeRaster(s3_ip_gg, "output/habitat_pred/s3_ip_gg.tif", overwrite = TRUE)

s3_ip_gg <- raster(x = "output/habitat_pred/s3_ip_gg.tif")

plot(s3_ip_gg)

registerDoMC(cores = 20)

s3_fp_gg <- foreach(i = seq_len(length(s3_fv_gg))) %dopar% {
  predict(s3_fv_gg[[i]],
          s3_mod_gg,
          type = "response",
          n.trees = s3_mod_gg$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s3_fp_gg-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s3_fp_gg <- stack(s3_fp_gg)

save(s3_fp_gg, file = "output/habitat_pred/s3_fp_gg")
writeRaster(s3_fp_gg, filename = "output/habitat_pred/s3_fp_gg.tif")

3.1.3 Scenario 4

s4_mod_gg <- gbm.step(data = s4_md_gg, gbm.x = 2:8, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s4_mod_gg)

# s4_ip_gg <- predict(s4_iv_gg, s4_mod_gg, type = "response", n.trees = s4_mod_gg$gbm.call$best.trees,  filename = "output/habitat_pred/4_fp_gg-000.tif")

# writeRaster(s4_ip_gg, "output/habitat_pred/s4_ip_gg.tif", overwrite = TRUE)

s4_ip_gg <- raster(x = "output/habitat_pred/s4_ip_gg.tif")

plot(s4_ip_gg)

registerDoMC(cores = 20)

s4_fp_gg <- foreach(i = seq_len(length(s4_fv_gg))) %dopar% {
  predict(s4_fv_gg[[i]],
          s4_mod_gg,
          type = "response",
          n.trees = s4_mod_gg$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s4_fp_gg-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s4_fp_gg <- stack(s4_fp_gg)

save(s4_fp_gg, file = "output/habitat_pred/s4_fp_gg")
writeRaster(s4_fp_gg, filename = "output/habitat_pred/s4_fp_gg.tif")

3.2 Leadbeater’s Possum

3.2.1 Scenario 1

s1_mod_lb <- gbm.step(data = s1_md_lb, gbm.x = 2:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s1_mod_lb)

# s1_ip_lb <- predict(s1_iv_lb, s1_mod_lb, type = "response", n.trees = s1_mod_lb$gbm.call$best.trees, filename = "output/habitat_pred/s1_fp_lb-000.tif")

# writeRaster(s1_ip_lb, "output/habitat_pred/s1_ip_lb.tif", overwrite = TRUE)

s1_ip_lb <- raster(x = "output/habitat_pred/s1_ip_lb.tif")

plot(s1_ip_lb)

registerDoMC(cores = 20)

s1_fp_lb <- foreach(i = seq_len(length(s1_fv_lb))) %dopar% {
  predict(s1_fv_lb[[i]],
          s1_mod_lb,
          type = "response",
          n.trees = s1_mod_lb$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s1_fp_lb-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s1_fp_lb <- stack(s1_fp_lb)

save(s1_fp_lb, file = "output/habitat_pred/s1_fp_lb")
writeRaster(s1_fp_lb, filename = "output/habitat_pred/s1_fp_lb.tif")

3.2.2 Scenario 3

s3_mod_lb <- gbm.step(data = s3_md_lb, gbm.x = 2:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s3_mod_lb)

# s3_ip_lb <- predict(s3_iv_lb, s3_mod_lb, type = "response", n.trees = s3_mod_lb$gbm.call$best.trees, filename = "output/habitat_pred/s3_fp_lb-000.tif")

# writeRaster(s3_ip_lb, "output/habitat_pred/s3_ip_lb.tif", overwrite = TRUE)

s3_ip_lb <- raster(x = "output/habitat_pred/s3_ip_lb.tif")

plot(s3_ip_lb)

registerDoMC(cores = 20)

s3_fp_lb <- foreach(i = seq_len(length(s3_fv_lb))) %dopar% {
  predict(s3_fv_lb[[i]],
          s3_mod_lb,
          type = "response",
          n.trees = s3_mod_lb$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s3_fp_lb-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s3_fp_lb <- stack(s3_fp_lb)

save(s3_fp_lb, file = "output/habitat_pred/s3_fp_lb")
writeRaster(s3_fp_lb, filename = "output/habitat_pred/s3_fp_lb.tif")

3.2.3 Scenario 4

s4_mod_lb <- gbm.step(data = s4_md_lb, gbm.x = 2:7, gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, step.size = 1, bag.fraction = 0.5, prev.stratify = FALSE, verbose = FALSE, max.trees = 1000)

summary(s4_mod_lb)

# s4_ip_lb <- predict(s4_iv_lb, s4_mod_lb, type = "response", n.trees = s4_mod_lb$gbm.call$best.trees, filename = "output/habitat_pred/s4_fp_lb-000.tif")

# writeRaster(s4_ip_lb, "output/habitat_pred/s4_ip_lb.tif", overwrite = TRUE)

s4_ip_lb <- raster(x = "output/habitat_pred/s4_ip_lb.tif")

plot(s4_ip_lb)

registerDoMC(cores = 20)

s4_fp_lb <- foreach(i = seq_len(length(s4_fv_lb))) %dopar% {
  predict(s4_fv_lb[[i]],
          s4_mod_lb,
          type = "response",
          n.trees = s4_mod_lb$gbm.call$best.trees,
          filename = paste0("output/habitat_pred/s4_fp_lb-", sprintf("%03d", i), ".tif"),
          overwrite = TRUE)
}

s4_fp_lb <- stack(s4_fp_lb)

save(s4_fp_lb, file = "output/habitat_pred/s4_fp_lb")
writeRaster(s4_fp_lb, filename = "output/habitat_pred/s4_fp_lb.tif")

4 Population viability analysis

4.1 Greater Glider

gg_trans_mat <- matrix(c(0.00,0.00,0.25,
                         0.50,0.00,0.00,
                         0.00,0.85,0.85),
                       nrow = 3, ncol = 3, byrow = TRUE) # based on Possingham et al 1994
colnames(gg_trans_mat) <- rownames(gg_trans_mat) <- c('Newborn','Juvenile','Adult')

gg_stable_states <- abs( eigen(gg_trans_mat)$vectors[,1] / base::sum(eigen(gg_trans_mat)$vectors[,1]) ) 

4.1.1 Scenario 1

s1_hs_gg <- stack(s1_ip_gg, s1_fp_gg)

for(i in 1:51){
  s1_hs_gg[[i]][is.na(s1_hs_gg[[i]][])] <- 0
}

s1_hs_gg <- mask(s1_hs_gg, mask = ch_mask)

s1_hab_k_gg <- calc(s1_hs_gg[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s1_hab_k_gg) <- "carryingCapacity"

s1_gg_popN <- stack(replicate(ncol(gg_trans_mat), s1_hab_k_gg))

s1_gg_popN <- s1_gg_popN*gg_stable_states

s1_gg_idx <- which(!is.na(s1_hs_gg[[1]][]) & s1_hs_gg[[1]][] < 0.95)

s1_gg_pop <- s1_gg_popN
s1_gg_pop[!is.na(s1_gg_pop)] <- 0
s1_gg_pop[[1]][sample(s1_gg_idx, 10000)] <- 1
s1_gg_pop[[2]][sample(s1_gg_idx, 10000)] <- 1
s1_gg_pop[[3]][sample(s1_gg_idx, 10000)] <- 1

s1_gg_pop <- s1_gg_pop*ch_mask

names(s1_gg_pop) <- colnames(gg_trans_mat)

s1_gg_TotpopN <- sum(cellStats(s1_gg_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s1_gg_init_pop_size <- sum(s1_gg_pop)

s1_gg_landscape <- landscape(population = s1_gg_pop,
                          suitability = s1_hs_gg,
                          carrying_capacity = s1_hab_k_gg)

s1_gg_pop_dynamics <- population_dynamics(change = growth(transition_matrix = gg_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s1_gg_results_1_50 <- simulation(landscape = s1_gg_landscape,
                         population_dynamics = s1_gg_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s1_gg_results_1_50, file = "output/pva_results/s1_gg_results_1_50.Rds")

plot(s1_gg_results_1_50)

plot(s1_gg_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.1.2 Scenario 3

s3_hs_gg <- stack(s3_ip_gg, s3_fp_gg)

for(i in 1:51){
  s3_hs_gg[[i]][is.na(s3_hs_gg[[i]][])] <- 0
}

s3_hs_gg <- mask(s3_hs_gg, mask = ch_mask)

s3_hab_k_gg <- calc(s3_hs_gg[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s3_hab_k_gg) <- "carryingCapacity"

s3_gg_popN <- stack(replicate(ncol(gg_trans_mat), s3_hab_k_gg))

s3_gg_popN <- s3_gg_popN*gg_stable_states

s3_gg_idx <- which(!is.na(s3_hs_gg[[1]][]) & s3_hs_gg[[1]][] < 0.95)

s3_gg_pop <- s3_gg_popN
s3_gg_pop[!is.na(s3_gg_pop)] <- 0
s3_gg_pop[[1]][sample(s3_gg_idx, 10000)] <- 1
s3_gg_pop[[2]][sample(s3_gg_idx, 10000)] <- 1
s3_gg_pop[[3]][sample(s3_gg_idx, 10000)] <- 1

s3_gg_pop <- s3_gg_pop*ch_mask

names(s3_gg_pop) <- colnames(gg_trans_mat)

s3_gg_TotpopN <- sum(cellStats(s3_gg_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s3_gg_init_pop_size <- sum(s3_gg_pop)

s3_gg_landscape <- landscape(population = s3_gg_pop,
                          suitability = s3_hs_gg,
                          carrying_capacity = s3_hab_k_gg)

s3_gg_pop_dynamics <- population_dynamics(change = growth(transition_matrix = gg_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s3_gg_results_1_50 <- simulation(landscape = s3_gg_landscape,
                         population_dynamics = s3_gg_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s3_gg_results_1_50, file = "output/pva_results/s3_gg_results_1_50.Rds")

plot(s3_gg_results_1_50)

plot(s3_gg_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.1.3 Scenario 4

s4_hs_gg <- stack(s4_ip_gg, s4_fp_gg)

for(i in 1:51){
  s4_hs_gg[[i]][is.na(s4_hs_gg[[i]][])] <- 0
}

s4_hs_gg <- mask(s4_hs_gg, mask = ch_mask)

s4_hab_k_gg <- calc(s4_hs_gg[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s4_hab_k_gg) <- "carryingCapacity"

s4_gg_popN <- stack(replicate(ncol(gg_trans_mat), s4_hab_k_gg))

s4_gg_popN <- s4_gg_popN*gg_stable_states

s4_gg_idx <- which(!is.na(s4_hs_gg[[1]][]) & s4_hs_gg[[1]][] < 0.95)

s4_gg_pop <- s4_gg_popN
s4_gg_pop[!is.na(s4_gg_pop)] <- 0
s4_gg_pop[[1]][sample(s4_gg_idx, 10000)] <- 1
s4_gg_pop[[2]][sample(s4_gg_idx, 10000)] <- 1
s4_gg_pop[[3]][sample(s4_gg_idx, 10000)] <- 1

s4_gg_pop <- s4_gg_pop*ch_mask

names(s4_gg_pop) <- colnames(gg_trans_mat)

s4_gg_TotpopN <- sum(cellStats(s4_gg_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s4_gg_init_pop_size <- sum(s4_gg_pop)

s4_gg_landscape <- landscape(population = s4_gg_pop,
                          suitability = s4_hs_gg,
                          carrying_capacity = s4_hab_k_gg)

s4_gg_pop_dynamics <- population_dynamics(change = growth(transition_matrix = gg_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s4_gg_results_1_50 <- simulation(landscape = s4_gg_landscape,
                         population_dynamics = s4_gg_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s4_gg_results_1_50, file = "output/pva_results/s4_gg_results_1_50.Rds")

plot(s4_gg_results_1_50)

plot(s4_gg_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.2 Leadbeater’s Possum

lb_trans_mat <- matrix(c(0.00, 0.00, 1.10,
                         0.59, 0.00, 0.00,
                         0.00, 0.59, 0.79),
                       nrow = 3, ncol = 3, byrow = TRUE) 
colnames(lb_trans_mat) <- rownames(lb_trans_mat) <- c('Newborn','Juvenile','Adult')

lb_stable_states <- abs( eigen(lb_trans_mat)$vectors[,1] / base::sum(eigen(lb_trans_mat)$vectors[,1]) ) 

4.2.1 Scenario 1

s1_hs_lb <- stack(s1_ip_lb, s1_fp_lb)

for(i in 1:51){
  s1_hs_lb[[i]][is.na(s1_hs_lb[[i]][])] <- 0
}

s1_hs_lb <- mask(s1_hs_lb, mask = ch_mask)

s1_hab_k_lb <- calc(s1_hs_lb[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s1_hab_k_lb) <- "carryingCapacity"

s1_lb_popN <- stack(replicate(ncol(lb_trans_mat), s1_hab_k_lb))

s1_lb_popN <- s1_lb_popN*lb_stable_states

s1_lb_idx <- which(!is.na(s1_hs_lb[[1]][]) & s1_hs_lb[[1]][] < 0.95)

s1_lb_pop <- s1_lb_popN
s1_lb_pop[!is.na(s1_lb_pop)] <- 0
s1_lb_pop[[1]][sample(s1_lb_idx, 3000)] <- 1
s1_lb_pop[[2]][sample(s1_lb_idx, 3000)] <- 1
s1_lb_pop[[3]][sample(s1_lb_idx, 3000)] <- 1

s1_lb_pop <- s1_lb_pop*ch_mask

names(s1_lb_pop) <- colnames(lb_trans_mat)

s1_lb_TotpopN <- sum(cellStats(s1_lb_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s1_lb_init_pop_size <- sum(s1_lb_pop)

s1_lb_landscape <- landscape(population = s1_lb_pop,
                          suitability = s1_hs_lb,
                          carrying_capacity = s1_hab_k_lb)

s1_lb_pop_dynamics <- population_dynamics(change = growth(transition_matrix = lb_trans_mat,
                                                       global_stochasticity = 0.3),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 2000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s1_lb_results_1_50 <- simulation(landscape = s1_lb_landscape,
                         population_dynamics = s1_lb_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s1_lb_results_1_50, file = "output/pva_results/s1_lb_results_1_50.Rds")

plot(s1_lb_results_1_50)

plot(s1_lb_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.2.2 Scenario 3

s3_hs_lb <- stack(s3_ip_lb, s3_fp_lb)

for(i in 1:51){
  s3_hs_lb[[i]][is.na(s3_hs_lb[[i]][])] <- 0
}

s3_hs_lb <- mask(s3_hs_lb, mask = ch_mask)

s3_hab_k_lb <- calc(s3_hs_lb[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s3_hab_k_lb) <- "carryingCapacity"

s3_lb_popN <- stack(replicate(ncol(lb_trans_mat), s3_hab_k_lb))

s3_lb_popN <- s3_lb_popN*lb_stable_states

s3_lb_idx <- which(!is.na(s3_hs_lb[[1]][]) & s3_hs_lb[[1]][] < 0.95)

s3_lb_pop <- s3_lb_popN
s3_lb_pop[!is.na(s3_lb_pop)] <- 0
s3_lb_pop[[1]][sample(s3_lb_idx, 3000)] <- 1
s3_lb_pop[[2]][sample(s3_lb_idx, 3000)] <- 1
s3_lb_pop[[3]][sample(s3_lb_idx, 3000)] <- 1

s3_lb_pop <- s3_lb_pop*ch_mask

names(s3_lb_pop) <- colnames(lb_trans_mat)

s3_lb_TotpopN <- sum(cellStats(s3_lb_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s3_lb_init_pop_size <- sum(s3_lb_pop)

s3_lb_landscape <- landscape(population = s3_lb_pop,
                          suitability = s3_hs_lb,
                          carrying_capacity = s3_hab_k_lb)

s3_lb_pop_dynamics <- population_dynamics(change = growth(transition_matrix = lb_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s3_lb_results_1_50 <- simulation(landscape = s3_lb_landscape,
                         population_dynamics = s3_lb_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s3_lb_results_1_50, file = "output/pva_results/s3_lb_results_1_50.Rds")

plot(s3_lb_results_1_50)

plot(s3_lb_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))

4.2.3 Scenario 4

s4_hs_lb <- stack(s4_ip_lb, s4_fp_lb)

for(i in 1:51){
  s4_hs_lb[[i]][is.na(s4_hs_lb[[i]][])] <- 0
}

s4_hs_lb <- mask(s4_hs_lb, mask = ch_mask)

s4_hab_k_lb <- calc(s4_hs_lb[[1]], fun = function(x){if(is.na(x)) x else rbinom(prob = x, size = 3, n = 1)})

names(s4_hab_k_lb) <- "carryingCapacity"

s4_lb_popN <- stack(replicate(ncol(lb_trans_mat), s4_hab_k_lb))

s4_lb_popN <- s4_lb_popN*lb_stable_states

s4_lb_idx <- which(!is.na(s4_hs_lb[[1]][]) & s4_hs_lb[[1]][] < 0.95)

s4_lb_pop <- s4_lb_popN
s4_lb_pop[!is.na(s4_lb_pop)] <- 0
s4_lb_pop[[1]][sample(s4_lb_idx, 3000)] <- 1
s4_lb_pop[[2]][sample(s4_lb_idx, 3000)] <- 1
s4_lb_pop[[3]][sample(s4_lb_idx, 3000)] <- 1

s4_lb_pop <- s4_lb_pop*ch_mask

names(s4_lb_pop) <- colnames(lb_trans_mat)

s4_lb_TotpopN <- sum(cellStats(s4_lb_pop, 'sum', na.rm = TRUE)) # Get total population size to check sensible
s4_lb_init_pop_size <- sum(s4_lb_pop)

s4_lb_landscape <- landscape(population = s4_lb_pop,
                          suitability = s4_hs_lb,
                          carrying_capacity = s4_hab_k_lb)

s4_lb_pop_dynamics <- population_dynamics(change = growth(transition_matrix = lb_trans_mat,
                                                       global_stochasticity = 0.1),
                                       dispersal = fast_dispersal(
                                         dispersal_kernel = exponential_dispersal_kernel(
                                           distance_decay = 8000)),
                                       modification = NULL)

plan(sequential, workers = 1)
s4_lb_results_1_50 <- simulation(landscape = s4_lb_landscape,
                         population_dynamics = s4_lb_pop_dynamics,
                         habitat_dynamics = NULL,
                         timesteps = 50,
                         replicates = 1)


saveRDS(object = s4_lb_results_1_50, file = "output/pva_results/s4_lb_results_1_50.Rds")

plot(s4_lb_results_1_50)

plot(s4_lb_results_1_50[1], type = "raster", stages = 0, timesteps = c(1:5))